Práctica de laboratorio 1 Regresión lineal simple y multiple sobre el dataset California.

En este caso usaremos el dataset California, incluido en los datasets de la Herramienta KEEL. Para ello, lo primero que tendremos que hacer es leer el dataset de una manera un tanto especial.

california <- read.csv("datasets/california.dat", comment.char="@")
#Asignación manual
names(california) <- c("Longitude", "Latitude", "HousingMedianAge",
"TotalRooms", "TotalBedrooms", "Population", "Households", "MedianIncome", "MedianHouseValue")
head(california)
##   Longitude Latitude HousingMedianAge TotalRooms TotalBedrooms Population
## 1   -118.23    33.80               26        239           135        165
## 2   -122.46    37.71               39       2076           482       1738
## 3   -122.06    37.94               19       4005           972       1896
## 4   -122.87    38.68               32       4073           718       2053
## 5   -122.47    37.66               18       4172           806       3226
## 6   -122.50    37.77               52       2299           441       1252
##   Households MedianIncome MedianHouseValue
## 1        112       1.3333           187500
## 2        445       3.1958           232100
## 3        893       2.5268           235700
## 4        629       3.7352           228000
## 5        790       5.7535           297900
## 6        415       5.0562           336700

Una vez nuestro dataset está creado, podemos añadirlo al entorno con attach por comodidad a la hora de referenciar a las variables, aunque para facilitar la comprensión del mismo esto no lo haremos.

Visualización.

Por medio de algunos gráficos, vamos a intentar buscar relaciones entre nuestras variables. Lo primero será usar la funcion pairs par ver que variables pueden resultarnos mas interesantes.

pairs(california)

En base al anterior gráfico, podemos crear otros gráficos, con la siguiente función de manera sencilla.

temp <- california
plotY <- function(x,y) 
{
  plot(temp[,y]~temp[,x], xlab=paste(names(temp)[x]," X",x,sep=""), ylab=names(temp)[y])
}

Acorde a los gráficos visto con la funcion pairs podemos encontrar cierta correlación entre variables. Vamos a dibujar sus gráficos.

plotY(5,4)

plotY(7,4)

Como vemos hay correlación lineal, pero si atendemos a los datos, aporta más bien poco ya que son cosas triviales, como a mas personas dentro de la casa mas habitaciones, algo totalmente trivial.

Vamos a analizar con más detalle las relaciones entre las variables independienes y la variable dependiente de nuestro problema.

plotY(3,9)

plotY(4,9)

plotY(5,9)

plotY(6,9)

plotY(7,9)

plotY(8,9)

Vamos a probar si realizando una transformación logarítmica podemos encontrar mejores relaciones.

logcalifornia <- log(california[3:9])
pairs(logcalifornia)

Volvemos a generar los gráficos de nuestros datos en función a la variable dependiente.

plot(logcalifornia$MedianHouseValue~logcalifornia$HousingMedianAge)

plot(logcalifornia$MedianHouseValue~logcalifornia$TotalRooms)

plot(logcalifornia$MedianHouseValue~logcalifornia$TotalBedrooms)

plot(logcalifornia$MedianHouseValue~logcalifornia$Population)

plot(logcalifornia$MedianHouseValue~logcalifornia$Households)

plot(logcalifornia$MedianHouseValue~logcalifornia$MedianIncome)

Los datos no son muy reveladores, pero parece que tenemso una buena candidata para el ajuste lineal simple. Esta candidata es el MedianInCome.

Obtención de modelos lineales simples

Vamos a obtener el modelo para MedianInCome.

fit1=lm(california$MedianHouseValue~california$MedianIncome)
fit1
## 
## Call:
## lm(formula = california$MedianHouseValue ~ california$MedianIncome)
## 
## Coefficients:
##             (Intercept)  california$MedianIncome  
##                   45083                    41794

Ahora que hemos construido el modelo simple, vamos a obtener los resultados.

summary(fit1)
## 
## Call:
## lm(formula = california$MedianHouseValue ~ california$MedianIncome)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -540700  -55951  -16978   36979  434025 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              45083.4     1322.9   34.08   <2e-16 ***
## california$MedianIncome  41794.2      306.8  136.22   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 83740 on 20637 degrees of freedom
## Multiple R-squared:  0.4734, Adjusted R-squared:  0.4734 
## F-statistic: 1.856e+04 on 1 and 20637 DF,  p-value: < 2.2e-16
plot(california$MedianHouseValue~california$MedianIncome)
abline(fit1,col="red") 

confint(fit1)
##                            2.5 %   97.5 %
## (Intercept)             42490.34 47676.46
## california$MedianIncome 41192.79 42395.56

En la anterior salida, cabe comentar el valor idéntico entre Adjusted R-squared y R-squared, debido a que solo estamos usando un predictor para nuestro modelo, si usaramos varios, deberíamos decantarnos por el Adjusted R-squared, ya que R-squared no se ajustaría a la realidad del problema.

Vamos a probar a crear ahora el modelo con nuestra versión con transformación logarítmica, para ver los resultados.

fitlog=lm(logcalifornia$MedianHouseValue~logcalifornia$MedianIncome)
summary(fitlog)
## 
## Call:
## lm(formula = logcalifornia$MedianHouseValue ~ logcalifornia$MedianIncome)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.62308 -0.27357 -0.00515  0.26005  2.61591 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                11.071215   0.008278  1337.4   <2e-16 ***
## logcalifornia$MedianIncome  0.814534   0.006222   130.9   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4207 on 20637 degrees of freedom
## Multiple R-squared:  0.4537, Adjusted R-squared:  0.4537 
## F-statistic: 1.714e+04 on 1 and 20637 DF,  p-value: < 2.2e-16
plot(logcalifornia$MedianHouseValue~logcalifornia$MedianIncome)
abline(fitlog,col="red") 

confint(fitlog)
##                                 2.5 %     97.5 %
## (Intercept)                11.0549896 11.0874403
## logcalifornia$MedianIncome  0.8023391  0.8267289

Los datos constatan que nuestra transformación de los datos a logarítmica no ha sido acerdata. Por ello, descartaremos esta medida.

Por último, vamos a calcular manualmente el error residual, para prácticar en el acceso de los datos del modelo y comprobar si ofrece resultados similares a los de la función summary().

names(fitlog)
##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "xlevels"       "call"          "terms"         "model"
sqrt(sum(fitlog$residuals^2)/length(fitlog$residuals))
## [1] 0.420654

Antes de entrar a predecir con nuestro modelo casos nuevos, vamos a probar con las demás variables ya que es posible que alguna de ellas aunque en el proceso de EDA pareciera que no eran buenas si que lo sean.

fit2<-lm(california$MedianHouseValue~california$HousingMedianAge)
fit3<-lm(california$MedianHouseValue~california$TotalRooms)
fit4<-lm(california$MedianHouseValue~california$TotalBedrooms)
fit5<-lm(california$MedianHouseValue~california$Population)
fit6<-lm(california$MedianHouseValue~california$Households)

summary(fit2)
## 
## Call:
## lm(formula = california$MedianHouseValue ~ california$HousingMedianAge)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -214479  -85039  -25833   58351  318941 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 179123.57    1985.54   90.21   <2e-16 ***
## california$HousingMedianAge    968.36      63.47   15.26   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 114800 on 20637 degrees of freedom
## Multiple R-squared:  0.01115,    Adjusted R-squared:  0.01111 
## F-statistic: 232.8 on 1 and 20637 DF,  p-value: < 2.2e-16
plot(california$MedianHouseValue~california$HousingMedianAge)
abline(fit2,col="red") 

confint(fit2)   
##                                   2.5 %     97.5 %
## (Intercept)                 175231.7612 183015.378
## california$HousingMedianAge    843.9577   1092.769
summary(fit3)
## 
## Call:
## lm(formula = california$MedianHouseValue ~ california$TotalRooms)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -313528  -86440  -26697   55726  311793 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           1.882e+05  1.248e+03  150.71   <2e-16 ***
## california$TotalRooms 7.098e+00  3.649e-01   19.45   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 114400 on 20637 degrees of freedom
## Multiple R-squared:  0.018,  Adjusted R-squared:  0.01796 
## F-statistic: 378.4 on 1 and 20637 DF,  p-value: < 2.2e-16
plot(california$MedianHouseValue~california$TotalRooms)
abline(fit3,col="red") 

confint(fit3)   
##                              2.5 %       97.5 %
## (Intercept)           1.857040e+05 1.905979e+05
## california$TotalRooms 6.382378e+00 7.812796e+00
summary(fit4)
## 
## Call:
## lm(formula = california$MedianHouseValue ~ california$TotalBedrooms)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -214300  -87389  -27788   57332  300592 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              1.994e+05  1.301e+03 153.239  < 2e-16 ***
## california$TotalBedrooms 1.387e+01  1.905e+00   7.284 3.35e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 115300 on 20637 degrees of freedom
## Multiple R-squared:  0.002564,   Adjusted R-squared:  0.002516 
## F-statistic: 53.06 on 1 and 20637 DF,  p-value: 3.354e-13
plot(california$MedianHouseValue~california$TotalBedrooms)
abline(fit4,col="red") 

confint(fit4)   
##                                 2.5 %       97.5 %
## (Intercept)              196844.77285 201945.68640
## california$TotalBedrooms     10.14098     17.60791
summary(fit5)
## 
## Call:
## lm(formula = california$MedianHouseValue ~ california$Population)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -195391  -86889  -26883   58179  308217 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            2.104e+05  1.291e+03 163.006  < 2e-16 ***
## california$Population -2.510e+00  7.091e-01  -3.539 0.000402 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 115400 on 20637 degrees of freedom
## Multiple R-squared:  0.0006067,  Adjusted R-squared:  0.0005582 
## F-statistic: 12.53 on 1 and 20637 DF,  p-value: 0.0004019
plot(california$MedianHouseValue~california$Population)
abline(fit5,col="red") 

confint(fit5)   
##                               2.5 %        97.5 %
## (Intercept)           207904.795019 212965.570709
## california$Population     -3.899878     -1.119958
summary(fit6)
## 
## Call:
## lm(formula = california$MedianHouseValue ~ california$Households)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -224837  -86861  -27949   56985  303059 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           1.969e+05  1.319e+03 149.313   <2e-16 ***
## california$Households 1.989e+01  2.097e+00   9.487   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 115200 on 20637 degrees of freedom
## Multiple R-squared:  0.004342,   Adjusted R-squared:  0.004294 
## F-statistic:    90 on 1 and 20637 DF,  p-value: < 2.2e-16
plot(california$MedianHouseValue~california$Households)
abline(fit6,col="red") 

confint(fit6)   
##                              2.5 %       97.5 %
## (Intercept)           194336.59833 199506.69451
## california$Households     15.78218     24.00202

Obtención de modelos lineales múltiples

Como hemos podido constatar los modelos simples no son muy reveladores, y la transformacion logaritmica aunque mejora en algunos casos, hace que nuestra mejor variable para el proceso de regresión empeore mas de un 2%. Basandonos en que un buen modelo de regresión deberá incluir a esta variable, vamos a generar modelos lineales múltiples con el fin de mejorar los valores obtenidos.

fit1<-lm(california$MedianHouseValue~california$MedianIncome+california$TotalRooms)
summary(fit1)
## 
## Call:
## lm(formula = california$MedianHouseValue ~ california$MedianIncome + 
##     california$TotalRooms)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -541281  -55947  -17009   36994  433866 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             45289.3575  1406.8203   32.19   <2e-16 ***
## california$MedianIncome 41820.8713   313.0262  133.60   <2e-16 ***
## california$TotalRooms      -0.1173     0.2726   -0.43    0.667    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 83740 on 20636 degrees of freedom
## Multiple R-squared:  0.4735, Adjusted R-squared:  0.4734 
## F-statistic:  9278 on 2 and 20636 DF,  p-value: < 2.2e-16
fit2<-lm(california$MedianHouseValue~california$MedianIncome+california$HousingMedianAge)
summary(fit2)
## 
## Call:
## lm(formula = california$MedianHouseValue ~ california$MedianIncome + 
##     california$HousingMedianAge)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -596759  -53839  -14999   36715  446731 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 -10198.86    1915.54  -5.324 1.02e-07 ***
## california$MedianIncome      43170.03     298.37 144.686  < 2e-16 ***
## california$HousingMedianAge   1744.30      45.04  38.728  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 80850 on 20636 degrees of freedom
## Multiple R-squared:  0.5091, Adjusted R-squared:  0.5091 
## F-statistic: 1.07e+04 on 2 and 20636 DF,  p-value: < 2.2e-16

En este caso, como hemos comentado anteriormente, debemos fijarnos en el valor Adjusted R-Squared ya que tenemos dos valores en nuestro modelo. Hemos conseguido mejorar bastante respecto al anterior por lo que parece que el modelo está empezando a comportarse mejor.

fit3<-lm(california$MedianHouseValue~california$MedianIncome+california$TotalBedrooms)
summary(fit3)
## 
## Call:
## lm(formula = california$MedianHouseValue ~ california$MedianIncome + 
##     california$TotalBedrooms)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -533222  -55959  -16543   36370  441499 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              36702.68    1518.13   24.18   <2e-16 ***
## california$MedianIncome  41821.50     305.91  136.71   <2e-16 ***
## california$TotalBedrooms    15.38       1.38   11.15   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 83490 on 20636 degrees of freedom
## Multiple R-squared:  0.4766, Adjusted R-squared:  0.4765 
## F-statistic:  9395 on 2 and 20636 DF,  p-value: < 2.2e-16

Mejora algo respecto a la inclusión de solo el valor de Median Income, por lo que parece que es buena candidata para un modelo con las tres. Vamos a seguir probando.

fit3<-lm(california$MedianHouseValue~california$MedianIncome+california$Population)
summary(fit3)
## 
## Call:
## lm(formula = california$MedianHouseValue ~ california$MedianIncome + 
##     california$Population)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -544701  -55908  -16877   37610  437659 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             49116.7020  1508.8881  32.552  < 2e-16 ***
## california$MedianIncome 41802.4722   306.5991 136.342  < 2e-16 ***
## california$Population      -2.8521     0.5144  -5.545 2.98e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 83680 on 20636 degrees of freedom
## Multiple R-squared:  0.4742, Adjusted R-squared:  0.4742 
## F-statistic:  9307 on 2 and 20636 DF,  p-value: < 2.2e-16
fit4<-lm(california$MedianHouseValue~california$MedianIncome+california$Households)
summary(fit4)
## 
## Call:
## lm(formula = california$MedianHouseValue ~ california$MedianIncome + 
##     california$Households)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -532066  -55865  -16574   36192  441915 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             36683.31    1514.20   24.23   <2e-16 ***
## california$MedianIncome 41748.83     305.91  136.48   <2e-16 ***
## california$Households      17.17       1.52   11.29   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 83480 on 20636 degrees of freedom
## Multiple R-squared:  0.4767, Adjusted R-squared:  0.4766 
## F-statistic:  9398 on 2 and 20636 DF,  p-value: < 2.2e-16

Tenemos por tanto dos claras candidatas para nuestro modelo, estás son: MedianInCome y HousingMedianAge, por lo tanto, podremos construir modelos con mas variables añadiendo algunas a estas. Para ver cuales añadir y cuales quitar, podemos crear el modelo con todas y ver cuales podemos descartar según su valor de p-value.

fitfinal<-lm(MedianHouseValue~.,data = california)
summary(fitfinal)
## 
## Call:
## lm(formula = MedianHouseValue ~ ., data = california)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -563016  -43593  -11324   30320  804281 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -3.594e+06  6.254e+04 -57.462  < 2e-16 ***
## Longitude        -4.282e+04  7.130e+02 -60.056  < 2e-16 ***
## Latitude         -4.258e+04  6.733e+02 -63.239  < 2e-16 ***
## HousingMedianAge  1.156e+03  4.317e+01  26.782  < 2e-16 ***
## TotalRooms       -8.196e+00  7.882e-01 -10.397  < 2e-16 ***
## TotalBedrooms     1.134e+02  6.902e+00  16.434  < 2e-16 ***
## Population       -3.855e+01  1.079e+00 -35.728  < 2e-16 ***
## Households        4.843e+01  7.516e+00   6.444 1.19e-10 ***
## MedianIncome      4.025e+04  3.351e+02 120.127  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 69530 on 20630 degrees of freedom
## Multiple R-squared:  0.6371, Adjusted R-squared:  0.637 
## F-statistic:  4528 on 8 and 20630 DF,  p-value: < 2.2e-16

Parece que todas las variables son relevantes en funcion del p-value, el modelo ha mejorado bastante desde que creamos la priemera aproximación y quitar algún elemento hará que este probablemente empeore ya que todos son considerados relevantes. Vamos a probar con nuestro modelo de transformación logaritmica a ver si ha mejorado algo.

fitlog<-lm(MedianHouseValue~.,data = logcalifornia)
summary(fitlog)
## 
## Call:
## lm(formula = MedianHouseValue ~ ., data = logcalifornia)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.1225 -0.2266  0.0042  0.2218  3.3338 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      11.600077   0.041789  277.59   <2e-16 ***
## HousingMedianAge  0.146326   0.004873   30.03   <2e-16 ***
## TotalRooms       -0.748241   0.015770  -47.45   <2e-16 ***
## TotalBedrooms     0.684295   0.021892   31.26   <2e-16 ***
## Population       -0.358934   0.009884  -36.31   <2e-16 ***
## Households        0.455336   0.019867   22.92   <2e-16 ***
## MedianIncome      1.107409   0.008141  136.04   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3725 on 20632 degrees of freedom
## Multiple R-squared:  0.5717, Adjusted R-squared:  0.5715 
## F-statistic:  4589 on 6 and 20632 DF,  p-value: < 2.2e-16

Podemos ver que comparándolo con el anterior modelo, este no ofrece tan buenos resultados. Por último calcularemos el error en nuestro modelo.

Obtención de error en test de nuestro modelo

Para ascercarnos más a la realidad de un problema dado, crearemos dos particiones del dataset california, una muestra aleatoria y del 20% para test y un 80% para el training.

sample <- sample.int(n = nrow(california), size = floor(.75*nrow(california)), replace = F)
train <- california[sample, ]
test  <- california[-sample, ]

Ahora con nuestro conjunto de entrenamiento creamos el modelo:

trainmodel<-lm(MedianHouseValue~.,data = train)
summary(trainmodel)
## 
## Call:
## lm(formula = MedianHouseValue ~ ., data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -556163  -43311  -10816   30524  473558 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -3.600e+06  7.205e+04 -49.968  < 2e-16 ***
## Longitude        -4.302e+04  8.207e+02 -52.418  < 2e-16 ***
## Latitude         -4.293e+04  7.747e+02 -55.411  < 2e-16 ***
## HousingMedianAge  1.120e+03  4.963e+01  22.565  < 2e-16 ***
## TotalRooms       -6.723e+00  9.027e-01  -7.448 9.99e-14 ***
## TotalBedrooms     1.062e+02  8.139e+00  13.048  < 2e-16 ***
## Population       -4.433e+01  1.364e+00 -32.492  < 2e-16 ***
## Households        6.154e+01  9.107e+00   6.757 1.46e-11 ***
## MedianIncome      3.958e+04  3.812e+02 103.830  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 69310 on 15470 degrees of freedom
## Multiple R-squared:  0.6389, Adjusted R-squared:  0.6387 
## F-statistic:  3421 on 8 and 15470 DF,  p-value: < 2.2e-16

Por último realizamos predicciones sobre el conjunto de test, que no hemos usado para probar el modelo y calculamos el error de manera manual.

yprime=predict(trainmodel,test)
sqrt(sum(abs(test$MedianHouseValue-yprime)^2)/length(yprime))
## [1] 70391.05

Podemos ver como el error residual estándard no está muy por debajo que los del conjunto de training, y aunque no es un valor muy bueno, es bastante mejor que cuando comenzamos con nuestro proceso de regresión lineal simple.

Práctica de laboratorio 2: Knn y Validación Cruzada

Vamos a usar knn sobre nuestro conjunto de datos, para ello si no lo tenemos debemos instalar el paquete kknn.

install.packages("kknn")
## 
## The downloaded binary packages are in
##  /var/folders/hx/tzq4mbtj1pj4gvxnfzdmx14h0000gn/T//RtmpN8yVkx/downloaded_packages
require(kknn)
## Loading required package: kknn

Una vez instalado el paquete, crearemos el modelo sobre nuestro conjunto de training.

fitknn1<-kknn(MedianHouseValue~., train, test)

Una vez creado el modelo, podemos ver los valores si accedemos a fitknn1$fitted.values. Para ver manualmente el error podemos realizar lo siguiente:

plot(test$MedianHouseValue~test$MedianIncome)
points(test$MedianIncome, fitknn1$fitted.values, col="blue", pch=20)

Podemos ver como la prediccion se ajusta bastante bien sobre lo0s datos reales, por lo que graficamente hemos comprobado el funcionamiento del modelo. También podemos caluclar el valor del error, para ello:

yprime<-fitknn1$fitted.values
sqrt(sum((test$MedianHouseValue-yprime)^2)/length(yprime))
## [1] 62172

El resultado es algo mejor que el de nuestro modelo lineal múltiple. Pero quizá aún podamos mejorarlo realizando algún tipo de modificación sobre las variables independientes.

fitknn2<-kknn(MedianHouseValue~.-HousingMedianAge, train, test)
yprime<-fitknn2$fitted.values
sqrt(sum((test$MedianHouseValue-yprime)^2)/length(yprime))
## [1] 60745.32

Hemos eliminado nuestra variable HousingMedianAge, del modelo, y parece que funciona un poco mejor. Con esto queda constatado, que aunque tengamos un buen modelo entre dos variables entre si, no se debe descartar la opción de probar que pasaría con los resultados si eliminamos una de ellas.

Validación cruzada.

Para evitar que nuestro modelo sobre-entrene, vamnos a realizar un proceso de validacion cruzada.

setwd("~/Desktop/MASTER CIENCIA DE DATOS/introduccion-ciencia-datos/Regresion/datasets/5fold")
nombre <- "california"
run_lm_fold <- function(i, x, tt = "test") 
{
  file <- paste(x, "-5-", i, "tra.dat", sep="")
  x_tra <- read.csv(file, comment.char="@") 
  file <- paste(x, "-5-", i, "tst.dat", sep="")
  x_tst <- read.csv(file, comment.char="@")
  In <- length(names(x_tra)) - 1 
  names(x_tra)[1:In] <- paste ("X", 1:In, sep="") 
  names(x_tra)[In+1] <- "Y"
  names(x_tst)[1:In] <- paste ("X", 1:In, sep="")
  names(x_tst)[In+1] <- "Y"
  if (tt == "train") 
  {
    test <- x_tra
  } 
  else
  {
    test <- x_tst
  }
  fitMulti=lm(Y~.,x_tra)
  yprime=predict(fitMulti,test)
  sum(abs(test$Y-yprime)^2)/length(yprime)
}
lmMSEtrain<-mean(sapply(1:5,run_lm_fold,nombre,"train")) 
lmMSEtest<-mean(sapply(1:5,run_lm_fold,nombre,"test"))

lmMSEtrain
## [1] 4826189710
lmMSEtest
## [1] 4844365688

Ahora lo aplicaremos para el algoritmo KNN:

setwd("~/Desktop/MASTER CIENCIA DE DATOS/introduccion-ciencia-datos/Regresion/datasets/5fold")
nombre <- "california"
run_knn_fold <- function(i, x, tt = "test") 
{
  file <- paste(x, "-5-", i, "tra.dat", sep="")
  x_tra <- read.csv(file, comment.char="@") 
  file <- paste(x, "-5-", i, "tst.dat", sep="")
  x_tst <- read.csv(file, comment.char="@")
  In <- length(names(x_tra)) - 1 
  names(x_tra)[1:In] <- paste ("X", 1:In, sep="") 
  names(x_tra)[In+1] <- "Y"
  names(x_tst)[1:In] <- paste ("X", 1:In, sep="")
  names(x_tst)[In+1] <- "Y"
  if (tt == "train") 
  {
    test <- x_tra
  } 
  else
  {
    test <- x_tst
  }
  fitMulti=kknn(Y~.,x_tra,test)
  yprime=fitMulti$fitted.values
  sum(abs(test$Y-yprime)^2)/length(yprime)
}
knnMSEtrain<-mean(sapply(1:5,run_knn_fold,nombre,"train")) 
knnMSEtest<-mean(sapply(1:5,run_knn_fold,nombre,"test"))
knnMSEtrain
## [1] 1560868807
knnMSEtest
## [1] 3845914481